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3 . We study the phase diagram of strongly interacting matter in the framework of a non- 

local SU(2) chiral quark model which includes wave function renormalization and coupling 
^vq ' to the Polyakov loop. Both non-local interactions based on the frequently used exponential 

form factor, and on fits to the quark mass and renormalization functions obtained in lattice 

re- 
calculations are considered. Special attention is paid to the determination of the critical 

■ points, both in the chiral limit and at finite quark mass. In particular, we study the position 

' of the Critical End Point as well as the value of the associated critical exponents for different 
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model parameterizations. 

PACS numbers: 12.39. Ki, 11.30. Rd, 12.38.Mh 



I. INTRODUCTION 

At low temperatures and densities strongly interacting matter is believed to be in a phase in 

which chiral symmetry is broken and the quarks are confined. However, as the temperature (T) 
• i-^ . 

and/or the chemical potential (fi) increase some kind of transition to a chiral restored and/or 

5-H ' 



deconfined phase is expected to happen. The detailed understanding of this phenomenon has 
become an issue of great interest in recent years, both theoretically and experimentallyjlj]. From 
the theoretical point of view, even if a significant progress has been made in the development of ab 
initio calculations such as lattice QCD 2- Jj, these are not yet able to provide a full understanding 
of the QCD phase diagram, due to the well-known difficulties of dealing with small current quark 
masses and finite chemical potentials. Thus, it is important to develop effective models that 
show consistency with lattice results and can be extrapolated into regions not accessible by lattice 
calculation techniques. Among them, the local Nambu-Jona-Lasinio (NJL) has been widely used 
to describe the behavior of strongly interacting matter at finite temperature and density [5|. In 
recent years an extension of the NJL model has been proposed in which the coupling of the 
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quarks 
model 



to the Polyakov loop is included. This so-called Polyakov-Nambu-Jona-Lasinio (PNJL) 



6Ml2l] allows to study the chiral and deconfmement transitions in a common framework. 
As an improvement over local models, the study of the phase diagram of chiral quark models that 
include nonlocal interactions 
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161 ] . These theories can be viewed as 



131 ] has been undertaken 
nonlocal extensions of the NJL model. In fact, nonlocality arises naturally in the context of several 
successful approaches to low-energy quark dynamics as, for example, the instanton liquid model 
\v\ and the Schwinger-Dyson resummation techniques 



3 3, Q 



also indicate that quark interactions should act over a certain range in momentum space. In 
addition, several studies 
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■ 24( have shown that nonlocal chiral quark models provide a satisfactory 
description of hadron properties at zero temperature and density. The aim of the present work is 
to extend previous studies of the chiral and deconfmement transitions in the framework of non- 



local chiral models with coupling to the Polyakov loop 
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2514271] by considering more general quark 



291 ] . we will adopt as the basic ingredient a reliable description 



interactions. Following Refs. 
of the T = n = quark propagator as given from fundamental studies, such as lattice QCD. In 
this sense, it should be noticed that most of the finite T and/or fj, calculations performed so far in 
the context of non-local chiral quark models have used exponential regulators and neglected the 
wave function renormalization (WFR) in the quark propagator. Recent lattice QCD calculations 
sug gest , however, that the WFR can be of the order of 30 % (or even more) at zero momentum 



19, 



201 ] . Moreover, these calculations also show that the quark masses tend to their asymptotic 



values in a rather soft way. Thus, it is important to perform a detailed study of the impact of the 
incorporation of these features on the predictions for the phase diagram and associated quantities. 
The lagrangian we will use is the minimal extension which allows to incorporate the full momentum 
dependence of the quark propagator, through its mass and wave function renormalization. Using 
such a model we will investigate the phase diagrams corresponding to different parameterizations, 
including one based on fits to the quark mass and renormalization functions obtained in lattice 
calculations, both in the chiral limit and for finite quark mass. The position of the critical points 
as well as the value of the associated critical exponents will be also studied. 

This article is organized as follows. In Sec. II we provide a description of the model and its 
parameterizations. In Sec. Ill we present and discuss the results obtained in the chiral limit, while 
those corresponding to finite values of the quark mass are given and analyzed in Sec. IV. In Sec. 
V we present a summary of our main results and conclusions. Finally, we include two appendices. 
In App. A we provide some details concerning the derivation of the Landau expansion associated 
with our model in the chiral limit, while in App. B we describe the formalism used to determine 
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the position of the Critical End Point. 

II. THE MODEL AND ITS PARAMETERIZATIONS 

We consider a nonlocal SU(2) chiral quark model which includes quark couplings to the color 
gauge fields. The corresponding Euclidean effective action is given by 

S E = J d A x U(x) (-i^D^ + m)^) -^j a (x)j a (x)-jp(x)j P (x)\ + U ($[A(x)])\ , (1) 

where ip is the Nf = 2 fermion doublet ip = (ii,d) T , and rh = diag(rn U: m^) is the current quark 
mass matrix. In what follows we consider isospin symmetry, that is m = m u = m^. The fermion 
kinetic term includes a covariant derivative = <9„ — iA^ where A^ are color gauge fields, and 
the operator J^d^ in Euclidean space is defined as 7 • V + 74 Jp, with 74 = i^Q. The nonlocal 
currents j a (x) , jp (x) are given by 



j a (x) = / d A z g(z) if) (x + r a if) (x - 



j P (x) = J d*z f(z) H^+IJ ^ 

Here, T a = (1,^75-?) and u(x') d v(x) = u(x')d x v(x) — d x iu(x')v(x). The functions g(z) and 
f(z) in Eq.([2]), are nonlocal covariant form factors characterizing the corresponding interactions. 
The four standard quark currents j a (x) require the same g(z) form factor to guarantee chiral 
invariance. The term ^jp(x)jp(x) is self-invariant under chiral transformations. The scalar- 
isoscalar component of the j a (x) current will generate the momentum dependent quark mass in 
the quark propagator, while the "momentum" current, jp(x), will be responsible for a momentum 
dependent wave function renormalization of this propagator. For convenience, we take the same 
coupling parameter, Gs, for both interaction terms. Note, however, that the relative strength 
between both interaction terms will be controlled by the mass parameter n p introduced in Eq.([2]). 
In what follows it is convenient to Fourier transform g(z) and f(z) into momentum space. Note 
that Lorentz invariance implies that the Fourier transforms g(p) and f(p) can only be functions 
of p 2 . 

To proceed we perform a standard bosonization of the theory. Thus, we introduce the bosonic 
fields cri t 2{x) and ir a (x), and integrate out the quark fields. In what follows, we work within 
the mean-field approximation (MFA), in which these bosonic fields are replaced by their vacuum 
expectation values a\ t 2 and ir a = 0. Since we are interested in studying the characteristics of 
the chiral phase transition we have to extend the so obtained bosonized effective action to finite 



4 



temperature T and chemical potential fi. In the present work this is done by using the Matsubara 
formalism. Concerning the gluon fields we will assume that they provide a constant background 
color field A4 = iAq = ig 5^0 Ga\ a /2, where Ga are the SU(3) color gauge fields. Then the traced 
Polyakov loop, which is taken as order parameter of confinement, is given by $ = ±Tr exp(z0/T), 
where <j> = iAq. We will work in the so-called Polyakov gauge, in which the matrix (f> is given a 
diagonal representation <fi = fyXa+faAg. This leaves only two independent variables, 03 and 0g. At 
vanishing chemical potential, owing to the charge conjugation properties of the QCD Lagrangian, 
the mean field traced Polyakov loop is expected to be a real quantity. Since 03 and 08 have to 
be real-valued 



31 



3Q |. this condition implies 



331 ] . As in e.g. Refs. 



10. 
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0. In general, this need not be the case at finite 



351 ] we will assume that the potential hi is such that the 
condition 08 = is well satisfied for the range of values of fi and T investigated here. The mean 
field traced Polyakov loop is then given by <3? = <1>* = [1 + 2 cos(03/T)]/3. 



Within this framework the mean field thermodynamical potential il MFA results 



MFA 



■1 p,n 



(p c -) 2 + M 2 (p c - 



+ 



2G S 



+ W($,T). 



(3) 



Here, the shorthand notation f pn = J2 n J d 3 p/ (27r) 3 has been used, and M(p) and Z{p) are given 
by 



Af(p) = Z(p)[m + a 1 g(p)} 
Z(p) = [1 - a 2 /(p)]- 1 . 



(4) 



In addition, we have defined 



Pn,: 



(2n + l)vrT -ifi + . 



(5) 



where the quantities C are given by the relation = diag(0 r , fl) Namely, C = c 03 with 
c = 1, — 1, for r, g, b, respectively. 

To proceed we need to specify the explicit form of the Polyakov loop effective potential. Fol- 
lowing Ref. [l(| we consider 

«($,T) = 



- a{T) <5> 2 + b(T) ln(l - 6 3> 2 + 8 $ 3 - 3 $ 4 ) 



where the coefficients are parametrized as 
a(T) = a + a 1 (^) + a 2 



b(T) = 63 



(6) 



(7) 



and the values of To, cij and 63 are fitted to QCD lattice results. 
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f2 MFA turns out to be divergent and, thus, needs to be regularized. For this purpose we use the 
same prescription as in e.g. Ref. Namely, 

fiMFA = n MFA _ n free + Q free + Qq (g) 

where Q,f ree is obtained from Eq.© by setting o\ = 02 = and Ureg B is the regularized expression 
for the quark thermodynamical potential in the absence of fermion interactions, 

= ~ iT j (0 E K 1 + -P ) + in (l + exp [-S±£±i*=] )" 

"(9) 

with E p = \J ft 1 + m 2 . Finally, note that in Eq.([8]) we have included a constant f^o which is fixed 
by the condition that vanishes at T = p = 0. 

The mean field values and <3? at a given temperature or chemical potential, are obtained 
from a set of three coupled "gap" equations. This set of equations follows from the minimization 
of the regularized thermodynamical potential, that is 

anMFA anMFA anMFA 

_ L reg _ u ^ L reg _ u ^ L reg _ „ ,-.„.. 

dai 8a 2 d$ 1 ' 

Once the mean field values are obtained, the behavior of other relevant quantities as a function 
of temperature and chemical potential can be determined. Here, we will be particularly interested 
in the chiral quark condensate (qq) and the quark density p defined by 

flnMFA anMFA 

as well as their corresponding susceptibilities, i.e. the chiral susceptibility Xch and the quark 
number susceptibility Xq-, defined by 

d(qq) dp 

Xch = ~d^ 5 x « = Tp- (12) 

Finally, the specific heat Cy, is expressed as 

Cv = -T . (13) 

In order to fully specify the model under consideration we have to fix the model parameters as 



well as t 
Ref. 



re form factors g(q) and f(q) which characterize the non-local interactions. Here, following 



29l | we consider two different types of functional dependencies for these form factors. The 



first one corresponds to the often used exponential forms, 

g(q) = exp {-q 2 /^l) \ f(q) = exp (V/A?) . (14) 



Note that the range (in momentum space) of the nonlocality in each channel is determined by the 
parameters Ao and Ai , respectively. Fixing the current quark mass and chiral quark condensate at 
T = fx = to the reasonable values m = 5.7 MeV and (qq) 1 ^ = 240 MeV, the rest of the parameters 
are determined so as to reproduce the empirical values f n = 92.4 MeV and m n = 139 MeV, and 
Z(0) = 0.7 which is within the range of values suggested by recent lattice calculations [8, 10]. In 
what follows this choice of model parameters and form factors will be referred as parametrization 
Set B. The second type of form factor functional forms we consider is given by 

l + a z a> m f m {q) - m a z f z (q) 1 + a z 

l + a z f z {q) a m -ma z 1 + a z f z (q) 

where 

Uq) = [l+ (^/Ag) 372 ]" 1 , f z {q) = [1+ (^/A?)]" 572 - (16) 
As shown in Ref. 



29J, taking m = 2.37 MeV, a m = 309 MeV, a z = -0.3, A = 850 MeV and Ai = 
1400 MeV one can very well reproduce the momentum dependence of mass and renormalization 
functions obtained in lattice calculations as well as the physical values of m n and /„■. In what follows 
this choice of model parameters and form factors will be referred as parametrization Set C. Finally, 
in order to compare with previous studies where the wave function renormalization of the quark 
propagator has been ignored we consider a third parametrization (Set A). In such a case we take 
Z(p) = 1 (i.e. f(p) = 0) and an exponential parametrization for g(p). Such a model corresponds 



to the "Scheme II" discussed in Ref. [241] . from where we take the parameters corresponding to 
(qq) 1 / 3 = 240 MeV. The values of the model parameters for each of the chosen parameterizations 
are summarized in Table I. 



III. PHASE DIAGRAM IN THE CHIRAL LIMIT 



In order to investigate the details of the phase diagram of the non-local models under study it 
is convenient to consider first the chiral limit m = 0. In this limit, general considerations imply 
that for sufficiently small values of chemical potential the chiral restoration transition is of second 
order with the transition temperature T c decreasing as jjl increases. At a certain value of \i = \itcp 
the transition becomes of first order. The point in the T — fj, plane defined by (Ttcp, Htcp) 
corresponds to the so-called "tricritical point" (TCP). For values of /i > ^tcp the corresponding 
T c continues to decrease until it reaches zero. This marks the end of the critical line, /i c (0) being 
the corresponding critical chemical potential. 
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In the following we will concentrate on the second order transition region. In such a region, for 
a given chemical potential fx, the condensate (qq) goes to zero when the temperature T approaches 
from below the critical value T c (/i), above which (qq) = and the chiral symmetry is restored. 
Thus, for T ~ T c (p,) the thermodynamical potential admits an expansion in powers of the order 
parameter (in this case the quark condensate). As discussed in detail in App. A, in the chiral limit 
such expansion reads 

n^ A = Clfa, T, * c , a 2c ) + A c (qq) 2 + C c (qq) 4 + O {(qq) 6 ) , (17) 

where the explicit expressions of and the coefficients A c and C c are given in Eqs. (|A3P and (]Alip . 
respectively. Having established the Landau expansion in terms of the chiral condensate as single 
independent variable, we can now analyze the characteristics of the phase transition following the 
standard textbook methods. For C c > the system undergoes a second order phase transition at 
a critical temperature T c . For each value of /i, this critical temperature can be obtained by solving 
a set of coupled equations given by the condition A c = supplemented by Eqs.( |A7|) . The values of 
T c (n) so obtained define a second order transition curve in the (T, /i) plane. As already mentioned, 
such a curve is a decreasing function of fi which starts at the critical temperature corresponding to 
vanishing chemical potential T c (0) and ends up at the tricritical point. The position of TCP can 
be determined by imposing the additional condition C c = 0. Namely, to obtain it one has to solve 
the set of coupled equations given by A c = C c = together with Eqs. (|A7D . 

To analyze the critical line beyond the TCP it is convenient to take T as independent variable 
and consider /U C (T). For T < Ttcp, the transition turns out to be discontinuous (i.e. first order). 
In this case, for each value of T, there is a region of values of \x for which three different solutions of 
the full gap equations, Eqs. ljlOp . exist. Two of them correspond to minima of the grand potential 
and the third one to a maximum. In the chiral limit considered in the present section, one of the 
minima has o\ = 0, while in the other o\ takes a finite (in general, non-negligible) value. Then, 
[i c corresponds to the chemical potential at which the pressure associated with these two minima 
coincide. 

The phase diagrams corresponding to our three parameterizations are displayed in Figfj] while 
the position of the characteristic points are given in Table [XTJ In Figfj] the dotted line indicates the 
second order chiral transition line, the full line that of first order and the dashed lines correspond 
to the deconfinement transition (the lower and upper lines correspond to ^ = 0.3 and $ = 0.5, 
respectively). We see that as \jl increases there appears a region where the system remains in 
its confined phase (signalled by $ smaller than ~ 0.3) even though chiral symmetry has been 
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restored. This corresponds to the recently proposed quarkyonic phase [361 ]. We observe that the 
general shape of the three diagrams is very similar with values of the critical temperatures at \i = 
differing by less than 4 MeV. In the case of the critical chemical potential at T = the difference 
between the three sets is somewhat larger. Comparing the result of Set A with that of Set B we see 
that the inclusion of the wave function renormalization implies a decrease of about 10 MeV in the 
value of /i c (0). The use of the softer form factors involved in the lattice inspired parametrization 
Set C leads to a further decrease of ~ 10 MeV . The feature of the phase diagram that turns out 
to be most sensitive to the model parametrization is the position of the TCP. In fact, although 
the three values of Ttcp are in a range of about 10 MeV, the value of fixcp increases by about 
a factor 2 when the wave function renormalization is included (i.e. when one goes from Set A to 
Set B) and by an extra factor ~ 3/2 when the lattice inspired parametrization Set C is used (i.e. 
when one goes from Set B to Set C). 

As it is well known, in the region of second order phase transition the behavior of several relevant 
thermodynamical quantities in the vecinity of the phase transition is determined by the critical 
exponents. In the case of the chiral and quark number susceptibilities, Xch and Xq respectively, 
and the specific heat Cy they are usually defined by 

Xch = \h - h c \-™ ; Xg = \h-h c \-^ ; C v = \h-h c \~ a , (18) 

where \h — h c \ is the distance to the critical point in the (//, T) plane. Note that in the chiral limit 
only trajectories approaching the transition from the chirally broken phase are relevant. In the 
present case, given the Landau expansion obtained above, one expects to have the usual mean field 
exponents. For trajectories which are not asymptotically tangential to the critical line they are 

Ich = 1 ; l q = ; a = 0, (19) 

for all points except for the TCP where 

lch = l ■ J q = l/2 ; a = 1/2. (20) 

As a test of consistency we have determined them numerically by studying the asymptotic behavior 
of the corresponding quantities for our three set of parameters. As an example of a typical result 
of such studies we show in the left panel of Figj2] the behavior of Xch f° r Set C as we approach 
an arbitrary point in second order transition line, i.e. a point different from the TCP, at constant 
H (for definiteness we consider [i = 10 MeV). The right panel in Fig{2] displays the results of a 
equivalent study for fj, = fiTCP- The values of the critical exponents extracted from this type of 



9 



analysis for Set C are given in Table IIIIl Very similar results are found for Set A and Set B. We 
see that in all cases the numerically obtained values are in very good agreement with the mean 
field ones given above. 

We finish this section by clarifying the role played by Polyakov loop in enhancing the critical 
temperature at a given value of fx, at least in the region where the transition is of second order. For 
simplicity we consider the parametrization Set A where there is no wave function renormalization. 
The condition A c = implies 1/8G = S21 (see App. A). However, following similar steps as those 
described in App.B of Ref . jl^ . it is possible to show that for T, fj, « Aq one has 



'21 - ^21 



Q app 

Or, 




3 





"3$ - r 




rp2 2\ 


^arccos 


2 


r 


T -n 1 



(21) 



In fact, we have checked numerically that in the relevant region T < 210 MeV and li < 50 MeV, 
Eq. (l2ip is verified with an accuracy higher than 15% for $ < 0.3. Therefore, the condition A c = 
leads to 

r c (/*) ^ c w - , (22) 



yJl-$ (arccos 



where 



Tc (fJ>) provides a good approximation to the critical temperature corresponding to pure quark 
(pq) non-local model, i.e. the model with no coupling to the PL, for the exponential regulator 
considered in parametrization Set A (see Ref. [3] for details). Of course, in the presence of PL- 



quark interactions the value of <l? c in Eq. (|22|) has to be obtained by simultaneously solving the 
corresponding gap equation, i.e. the second equation in Eq. (|A7p in the present case. However, 
we clearly see that for any value of $ c < 1 we have T c {ll) > Tj" J (/i). For example, at fx = one 
typically has $ c ~ 0.2 which implies T c (0) /T^ 9 (0) ~ 1.66. Since for the parametrization Set A 
we have T^ 9 (0) = 126 MeV for the pure quark non-local model in the chiral limit, we see that 
the coupling to the PL is expected to raise this value up to T c ~ 209 MeV which is in very good 
agreement with the numerically found value listed in Table ITU As it is clear from Eqs. ()22)23j) a 
similar enhancement of the critical temperature can be obtained at (low) finite 11. On the other 
hand, it should be noticed that in order to apply the present type of analysis to relate the values of 
li c (T) predicted in models with and without PL one must have a common range of temperatures 
for which the transition is of second order. However, for the parametrizations considered here this 
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is not possible since they always lead to Ttcp > Tj" J (0). For example, from Table II we see that 
Set A leads to Ttcp — 204.8 MeV to be compared with the value Tf q (0) = 126 MeV quoted above. 

IV. PHASE DIAGRAM FOR FINITE QUARK MASS 

We start by analyzing the behavior of some mean field quantities as functions of T and \x. 
Since the results obtained for our three different parameterizations are qualitatively quite similar 
we only present explicitly those corresponding to the parametrization Set C. They are given in 
Figj3] where we plot o\ , 02 and 3> as functions of T for some representative values of the chemical 
potential. The left panel of Figj3] shows that at fj, = there is a certain value of T at which o\ 
drops rapidly signalling the existence of a chiral symmetry restoration crossover transition. At 
basically the same temperature the Polyakov loop increases which can be interpreted as the 
onset of the deconfinement transition. As fi increases there is a certain value of fj, = [icep above 
which the transition starts to be discontinuous. At this precise chemical potential the transition is 
of second order. This situation is illustrated in the central panel of FigJSJ The corresponding values 
{TcePi ^cep) define the position of the so-called "critical end point" (CEP) which, as explained 
in App. B, can be found by solving a system of equations formed by the gap equations Eqs. (|10p 
supplemented by two additional equations of the type of Eqs. (|B3IB4|) . As displayed in the right 
panel of FigG2 for fi > ^cep the transition becomes discontinuous, i.e. of first order. Finally, for 
chemical potentials above [i c (T = 0) ~ 310 MeV the system is in the chirally restored phase for all 
values of the temperature. It is important to note that although 02 appears to be rather constant 
in Fig J3j at higher values of T it does go to zero as expected. 

The different nature of the chiral transition in each of the three regions of Figj3] is even more 
clearly observed in the behavior of the corresponding response functions. In FigJ5]we display the 
specific heat Cy as well as the chiral and quark number susceptibilities, Xch an d x<?; as a function 
of the temperature for parametrization Set C and the three different values of [i used in Figj3l 
The dotted line corresponds to [i = 0. We observe that all the response functions show a rather 
broad peak of finite height at basically the same value of T. Such a value of T corresponds to 
the temperature at which the crossover transition occurs. The dashed-dotted line corresponds to 
I 1 = Hcep which indicates that all the response functions display a sharp and narrow divergent 
peak. Such a behavior signals the second order nature of the chiral transition at the CEP. Finally, 
the full line corresponds to [i > ^cep- In ah cases we observe a discontinuity in the response 
functions which indicates that the associated transition is of first order. 
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The phase diagrams corresponding to our three different parameterizations are given in Figj5j 
Here the dotted line represents the line of crossover chiral transition while the full line that of first 
order. The dashed lines are associated to the deconfinement transition (the lower and upper lines 
correspond to $ = 0.3 and <3? = 0.5, respectively). The position of the most relevant points in the 
phase diagrams are tabulated in Table HVl As in the chiral case, we observe that the main difference 
appears in the position of the point at which the first order transition line ends. Comparing the 
results of Set A and Set B we see that the main effect of the wave function renormalization term 
is to shift the location of the CEP towards lower values of T and higher values of \i. Concerning 
the lattice motivated parametrization Set C we observe that it leads to even lower values of Tqep 
and higher values hcep- Comparing with the results obtained in the chiral limit we note that 
the variation of Tqep between different sets is larger in this case. Moreover, considering each 
parametrization separately we observe that while the effect of introducing a finite quark mass on 
the values of the T c (0) and // c (0) is quite small (aprox. +5 MeV and +30 MeV, respectively), the 
position of the CEP is much more sensitive to the value of the current quark mass m, especially 
in the case of parametrization Set C. This is clearly seen in FigE] where we plot Tqep and hcep 
as a function of m. In all the cases we note a sharp decrease (increase) of Tqep {h-cep) for low 
values of m. In the case of the exponential parametrization Set A and Set B, at m ~ 4 MeV this 
variation tends to disappear and the position of the CEP remains rather stable up to the maximum 
value of m we have considered. On the other hand, for the lattice motivated parametrization the 
situation is somewhat different. In fact, the variation is rather large for basically all the values 
of m considered. In particular we see that, after the initial sharp decrease, at about 4 MeV the 
value of Tqep starts to increase in a rather pronounced way. It is interesting to note that the 



behavior of the position of t 
be analytically investigated 



le CEP as a function of m close to m = (i.e. close to the TCP) can 



37] . General arguments imply that 



ATcep = T CE p(m) - T TCP = -c m 2 / 5 + 0(m 4 / 5 ) 

Aucep = HCEp{m) - htcp = +d m 2 / 5 + 0(m 4 / 5 ), (24) 

where c and d are definite positive constants. To verify that our results do satisfy these relations 
we have numerically study in detail the variation of Tqep and Hcep as a function of m close to 
the chiral limit. Results for the parametrization Set C are shown in Fig J7] where the power-law 
behavior of both ATcep and Aucep is clear seen as a straight line in the corresponding log-log 
plot. Performing a linear fit we obtain that the slope of both straight lines is 0.40 ± 0.01 in very 
good agreement with the exponents in Eq. fpT 
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Finally, we consider the behavior of the response functions close to the CEP. As already men- 
tioned, the chiral phase transition at this point is of second order. Thus, a critical behavior with 
critical exponents defined as in Eq. (jl8[) is expected. Within the approximations used in this work 
these exponents should take the mean field values 7^ = 7 g = a — 2/3 when one approaches the 
CEP using trajectories which are not asymptotically parallel to the first order transition line. As 
mentioned in Sec. Ill, these exponents can be numerically obtained by analyzing the asymptotic 
behavior of the response functions close to the critical point. In Figj8]we show a typical result of 
this type of numerical study. There, we display a log-log plot of the specific heat Cy corresponding 
to the parametrization Set C as a function of the relative temperature departure \T — Tcep\/Tcep 
for trajectories that approach the CEP at constant ^ = ^cep both from below (i.e. T < Tqep) 
and from above (i.e. T > Tqep)- We observe that a single straight line behavior is obtained up 
to relative departure as large as 10 -3 after which non-linear effects start to show up. This result 
is particularly interesting since there have been claims [371. |38|| that in some cases there might be a 



"two straight lines" behavior, namely, that two different critical exponents are needed in order to 
describe the critical behavior of the Cy close to the CEP. In fact, this feature has been interpreted 
as an influence of the TCP critical properties on the CEP ones. As displayed in Figj8] no sign 
of this type of effect is found in our case. Results similar to those presented in this figure have 
been obtained for the chiral and quark number susceptibilities. The corresponding values of the 
critical exponents for different kind of trajectories are listed in Table [V] Note that in all cases 
the agreement with the mean field values of the exponents is very good. The same type of results 
has been found for the other two parameterizations, i.e. the exponential parameterizations Set A 
and Set B. It should be stressed that in order to obtain numerically the critical exponents with 
a good accuracy it is important to know the position of the CEP with very good precision. We 
have checked that even an error of 0.5 MeV in Tqep and/or \xcep translates into an important 
uncertainty in the critical exponents. In this sense, the method for the determination of the CEP 
discussed in App. B is of great help. 



V. SUMMARY AND CONCLUSIONS 



In this work we have studied the behavior of strongly interacting matter at finite temperature 
and chemical potential using a non-local chiral quark model which includes wave function renor- 
malization and coupling to the Polyakov loop. This type of model can be understood as a non-local 
extension of the local Polyakov-Nambu-Jona-Lasinio model, and represents a step towards a more 
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realistic modeling of the QCD interactions that could allow a simultaneous description of the decon- 
finement and chiral phase transition. The non-local interactions have been described by considering 
not only the frequently used exponential form factors, but also a parametrization based on fits to 
the quark mass and renormalization functions obtained in lattice calculations. In this framework, 
we have studied the corresponding phase diagrams and associated quantities, both in the chiral 
limit and at finite values of the current quark mass, paying particular attention to the accurate 
determination of the critical points. In fact, in both cases we have been able to obtain a set of 
coupled equations for the position of the corresponding critical point, i.e. the Tricritical Point 
(TCP) in the chiral limit and the Critical End Point (CEP) for finite quark mass. Our numerical 
results indicate that some of the features of the phase diagrams are not very much dependent on 
the different parameterizations we used. For example, for finite quark mass the critical tempera- 
tures at fx = are within the range 210 — 215 MeV, while the critical chemical potentials at T = 
are in the range of 298 — 322 MeV. On the other hand, the position of the critical point turns 
out to be very sensitive to both the parametrization and the value of the current quark mass m. 
Comparing the results corresponding to the exponential parametrization with and without quark 
wave function renormalization we find that the main effect of the presence of this term is to shift 
the location of the CEP towards lower values of T and higher values /x. Concerning the lattice 
motivated parametrization we observe that it leads to even lower values of Tq e p an d higher values 
HCEP- As for the dependence on m we have verified that for small values of m (i.e. close to the 
TCP) the position of the CEP displays in all cases a power-law behavior, as expected. For the 
exponential parameterizations at m ~ 4 MeV this initial variation tends to disappear and the posi- 
tion of the CEP remains rather stable up to the maximum value of m we have considered. On the 
other hand, for the lattice motivated parametrization the situation is somewhat different. In fact, 
the dependence on m is rather strong for basically all the values considered. In particular, after an 
initial sharp decrease, at m ~ 4 MeV the value of Tqep starts to increase in a rather pronounced 
way. Finally, we have analyzed numerically the critical behavior around the TCP and the CEP 
determining the critical exponents associated with the chiral and quark number susceptibilities 
as well as the heat capacity. In all cases, we find that the obtained exponents agree with their 
predicted mean field values to a rather good degree of accuracy. In particular, no influence of the 
TCP properties on the CEP critical exponents has been found. 
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APPENDIX A: Derivation of the Landau expansion 



To derive the Landau expansion Eq. (|17p we start by assuming that the chemical potential ji is 
such that, in the chiral limit, the chiral condensate vanishes at a critical temperature T c (n). Since 
in that situation the mean field value v\ also vanishes, for T ~ T c (fi) it is possible to perform a 



double expansion of ^l^ A , Eq.©, in powers of v\ and m. We obtain 



n^ A (fi,T,^,a 2 ,ax) = £l(fi,T,<f>,v 2 ) + 4 



reg 



1 

8G 



S 2 i{fi,T,$,a 2 ) 



<j\ + 2 S 4 2(/i,r,$,<72) A 



-8ma 1 [S n ^,T,^,a 2 )-S 32 (fx,T,^ ) a 2 ) vf) + 0(vf,m vim 2 of) , 

(Al) 



where 



S jk (^T,^,v 2 ) =J2 I 9 3 (Pn,p) 

„ J v,n 



c \2 



(A2) 



and 0,(fi, T, a 2 ) is the MFA thermodynamical potential in the chiral limit for vanishing o\. 
Namely, 

8T 



kI v 2 



(A3) 

■ 

c J P,n 

Using Eq. (|lip the corresponding expression for the chiral condensate can be readily obtained. In 
the chiral limit we get 



(qq) = 8 vi [S n (n,T,$,v 2 ) - S 32 ^,T,^,v 2 )vj] +0{v\). 
Inverting this equation and replacing in Eq, (|Aip we finally get 



(A4) 



n^ g A (fi,T^,a 2 ,{qq)) = fi(/x, T, v 2 ) + A(ji, T, v 2 ) {qq)' 



+ C(/i,r,$,a 2 ) (qq) 4 + 0((qq)' 



(A5) 



Here, the coefficients A and C are given by 



A(n,T,$,v 2 ) 
C(fi,T,$,v 2 ) 



1 



4 ( S 1 2 1 (T, M ,<I>,<7 2 ) 
S A2 {T,n,$,a 2 ) 



±-S 21 (T,n,$,v 2 ) 
S 32 {T,^,v 2 ) 



128 Sf, (T, n, v 2 ) 32 S* (T, M , <D, <r 2 ) 



8G 



5 2 i(T,/i,$,c7 2 ; 



• (A6) 



15 



n 



The expansion Eq.(jA5]) looks very similar to Eq.(15) of Ref. [151] where non-local models in the 
absence of wave function renormalization and coupling to the Polyakov loop were analyzed. In 
principle, following similar ideas, explicit equations for the second order transition line as well as 
the position of the TCP point might be determined. However, the fact that in the present case the 
coefficients A and C depend on the mean field values a 2 and $ introduces further complications. 
In fact, slightly below T c (for a given value of //) the non- vanishing value of the condensate induces 
departures of ($,£12) from the critical values (3> C )f2,c) obtained as solutions of the set of equations 



an(/i,r c (/i),$,<7 2 ) 



da 2 



dn(ji,T c (ji),$,a2) 



<9$ 



0. 



(A7) 



Those departures can be obtained using the gap equations resulting from Eq. (jA5[) . To quadratic 
order in the condensate, we obtain 



(7 2 



91M (dl 2(T2 n c ) - (dl a n c 



di 2(T2 n c (a*40- [dLfic ){d n A 



91M (dl 2a2 n c 



9LA 



(A8) 



where a compact notation has been used to denote the derivatives of Ci and A evaluated at T c . In 
this notation we have, for example, 



0fi(/i,r c Oi),$,ff 2 ) 



d<S> 



dA(fx,T c (n),^,a 2 ) 



9(72 



(A9) 



*c<T2,c da 2 
Using Eqs. (|A8p we can now obtain the leading corrections to A, C and Q induced by the non- 
vanishing value of the condensate. Replacing the corresponding results in Eq. (|A5p . and grouping 
in powers of (qq) we finally obtain Eq. (|17p . Namely, 



qMFA 
iL reg 



T, $ c , a 2)C ) + A c (qqf + C c (qq)* + O {(qqf) , 



(A10) 



where 



A c = i4(/x,r c (/x),$ c ,er 2)C ), 
C c = C(n,T c (n),^ c ,a 2t c) - 

d 2 a2a2 n c ) (ds,A c f + (dl h c ) (d a2 A c f - 2 (dl a h c ) (d U2 A c ) (d a2 A c 



(dl 



<£><J> J "C 



(dl n „n, 



SLA 



(All) 
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APPENDIX B: Formalism to determine the position of the CEP 



As discussed in Ref. [37], in the case in which the grand potential depends on one single 
variational parameter, a set of equations that allows to determine the position of the CEP can 
be obtained. This set is formed by the corresponding gap equation supplemented with the two 
equations that result from demanding that the second and third of grand potential with respect 
to the parameter also vanish. Physically, this corresponds to determining the values of (T, fj,) for 
which the grand potential around its minimum is as flat as possible. The purpose of this appendix 
is to give some details of the formalism needed to generalize this idea to the case in which the 
grand potential depends on more than one variational parameter. In order to keep the derivation 
as general as possible we will assume here that the grand potential depends on iV variational 
parameters £i, £2, •■> £,N- For convenience, in what follows, we will distinguish £1 from the rest (note 
that there is no loss of generality in this choice since the ordering of the parameters is completely 
arbitrary), and introduce the latin index j = 2, N . Then, in a compact notation, the set of gap 
equations reads 

f - ° 

The basic idea is now to consider the variational parameters £j as functions of £1 with the corre- 
sponding functional dependence determined by the solutions of the N — 1 gap equations Eqs. (|B2p , 
The total set of N + 2 equations needed to determine the CEP is then obtained by supplementing 
these N — 1 equations with those resulting from demanding that the first, second and third total 
derivatives of the grand potential with respect to £1 vanish. The first of these equations turns out 
to be, of course, the gap equation Eq. (|Bl|) . The other two result 



d 2 Q d 2 Q d 2 Q, 

~atT + 2 a? a? £?' + at at £7' ^'k = °> ( B3 ) 



d 3 n d 3 n , d 3 n , , 

f 3 M2ac. £j + 3 atacat, £fc 



Here, and in what follows, the sum over repeated latin indexes k,j, I = 2, .., iV is understood. Note 
that in obtaining these equation some terms have been dropped assuming that the gap equations 
Eqs. (IBl)B2j) are simultaneously satisfied. In Eqs. (IB3|B4|) . £j , etc., denote the derivatives of 
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the corresponding parameters with respect to £1. They can be conveniently expressed in terms of 
partial derivatives of the grand potential by solving the two set of linear equations resulting from 
taking the first and second total derivatives of both sides of the gap equations Eq. (|B2jl . We obtain 

d 2 n 



C'j = -(cr l ) jk (B5) 



= -(c- 1 ) ifc 



(B6) 



where 



Therefore, once the explicit form of the grand potential O in terms of the variational parameters 
is known, all the derivatives appearing in the set of equations Eqs. (|Bll - IB4p can be analytically 
determined. Then, the numerical solution of this set of equations allows to determine the values of 
(T, jj) corresponding to the CEP as well as the corresponding values of the variational parameters 

Turning to the model discussed in the present work, we note that for the parameterizations Set 
B and Set C we have three variational parameters. We chose to identify them as £1 = °i> £2 = a 2 
and £3 = <&. Then, the set of five equations needed to determine the CEP is formed by Eqs. (|10p 
supplemented by the two equations that result from performing the corresponding identifications 
in Eqs. (IB3l - IB7p , In the case of parametrization Set A the situation is somewhat simpler since 
there are only two variational parameters (£1 = o\ and £2 = 3>) and, hence, four equations are 
required to determine the CEP. 
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TABLE I: Set parameters and chiral condensates for T = \i = 0. 







Set A 


Set B 


Set C 


m c 


MeV 


5.78 


5.70 


2.37 


G s Kl 




20.650 


32.030 


20.818 


A 


MeV 


752.20 


814.42 


850.00 


Hp 


GeV 




4.180 


6.034 


Ai 


MeV 




1034.5 


1400.0 


0"1 


MeV 


424 


529 


442 


(72 






-0.43 


-0.43 


- < qq > 1/3 


MeV 


240 


240 


326 



TABLE II: Position of some characteristic points of the phase diagrams in the chiral limit. All values are in 
MeV. 





Set A 


Set B 


Set C 


T c (0) 


206.6 


205.9 


209.7 


\±TCP 


46.2 


86.1 


125.7 


Ttcp 


204.8 


199.2 


194.6 


Mc(0) 


297.6 


285.5 


268.2 
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TABLE III: Critical exponents in the chiral limit for Set C. 







Ich 


7? 


Q 


Point in 




1.00(1) 


0.00(1) 


0.00(1) 


2nd order 


Tt 


1.00(1) 


0.00(1) 


0.00(1) 


critical line 


MF exponent 


1 












1.00(1) 


0.51(1) 


0.50(1) 


TCP 


Tt 


1.00(1) 


0.51(1) 


0.50(1) 




MF exponent 


1 


1/2 


1/2 



TABLE IV: Position of some characteristic points of the phase diagrams for finite quark mass. All values 
are in MeV. 





Set A 


Set B 


Set C 


T c (0) 


210.0 


209.8 


214.5 


P-CEP 


132.5 


182.3 


234.8 


Tqep 


197.8 


181.6 


154.2 




321.5 


311.6 


298.1 



TABLE V: Critical exponents at the CEP for Set C. 





Ich 




n 




0.67(1) 


0.67(1) 


0.66(1) 


n <- 


0.66(1) 


0.66(1) 


0.67(1) 




0.67(1) 


0.67(1) 


0.66(2) 


Tl 


0.66(1) 


0.66(1) 


0.67(1) 


MF exponent 


2/3 


2/3 


2/3 
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FIG. 1: Phase diagrams in the chiral limit for the three paramcterizations considered. Set B and Set 
C include quark wave function rcnormalization while Set A does not. Set A and Set B correspond to 
exponential form factors while Set C to lattice motivated form factors. The dotted line corresponds to the 
second order chiral transition and the full line to that of first order one. The dashed lines correspond to the 
deconfmcmcnt transition (the lower and upper lines being for $ = 0.3 and $ = 0.5, respectively). 




FIG. 2: Dependence of the chiral susceptibility Xch as a function of T for constant fi in the vicinity of an 
arbitrary point (taken to correspond to \x = 10 MeV) in the 2nd order transition line (left panel) and the 
TCP (right panel) in the chiral limit for parametrization Set C. 
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FIG. 3: Mean fields oi, a-i and $ as functions of T at three representative values of chemical potentials for 
parametrization Set C. Note that the scale to the left corresponds to u\ while that to the right to ai and 
i>. Since 02 turns out to be negative we plot — a^. 
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FIG. 4: Behavior of the specific heat CV, the chiral susceptibility Xch and the quark number susceptibility 
X q as functions of T at three representative values of chemical potentials for parametrization Set C. 
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FIG. 5: Phase diagrams for the three parameterizations considered. Set B and Set C include quark wave 
function renormalization while Set A does not. Set A and Set B correspond to exponential form factors while 
Set C to lattice motivated form factors. The dotted line corresponds to the line of crossover chiral transition 
and the full line to that of first order chiral transition. The dashed lines correspond to the deconfinement 
transition (the lower and upper lines being for $ = 0.3 and <E> = 0.5, respectively). 
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FIG. 6: Position of the CEP in the (T,/j,) plane as a function of the current quark mass m for the three 
model parameterizations used in this work. The left panel displays the dependence of Tcep while the right 
panel that of \icep- Values corresponding to the chiral limit m = and to the "physical" current quark 
masses given in Table I for each parameter set are indicated by fat dot. 
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FIG. 7: Position of the CEP in the (T, /x) plane as a function of the current quark mass m close to the chiral 
limit for parametrization Set C. 




FIG. 8: Dependence of the specific heat Cy as a function of T for constant /i in the vicinity of the CEP 
for parametrization Set C. The left panel displays the dependence as one approaches the CEP from below 
while the right panel the one when the approach is done from above. 



